//=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|
// Classe GaussianDistribution:
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//    Usa a funcao do NR para gerar um numero randomico com distribuicao 
//  gaussiana de media zero e varianca 1.
//    Por default, a semente eh time(0)*PID.
//________________________________________________________
//A.Z. - 01/06 => Criacao
//       02/06 => Usa a biblioteca seed.h
//=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|
#ifndef GAUSS_H
#define GAUSS_H

#include <ctime>
#include <cmath>
#include "seed.h"
#include "meusTipos.h"
using namespace std;



class GaussianDistribution{
public:

  // O construtor garante que a semente seja NEGATIVA na primeira chamada
  GaussianDistribution( MTci& seed = mySeed ){ _seed = -abs(seed); }


  // Retorna um numero randomico
  inline double operator () (void){
    return _NRCPP_gasdev(_seed);
  }

private:
  inline double _NRCPP_gasdev(int &);
  inline double _NRCPP_ran1(int &);

  int _seed;
};





//------------------------------------------------------------------------------
// DESCRICAO:
//   Funcoes do NR C++
inline double GaussianDistribution::_NRCPP_gasdev(int &idum){
  static int iset=0;
  static double gset;
  double fac,rsq,v1,v2;

  if(idum<0) iset=0;
  if(iset==0){
    do{
      v1=2.0*_NRCPP_ran1(idum)-1.0;
      v2=2.0*_NRCPP_ran1(idum)-1.0;
      rsq=v1*v1+v2*v2;
    }while(rsq>=1.0 || rsq==0.0);
      fac=sqrt(-2.0*log(rsq)/rsq);
      gset=v1*fac;
      iset=1;
      return v2*fac;
  }else{
    iset=0;
    return gset;
  }
}
inline double GaussianDistribution::_NRCPP_ran1(int &idum){

  MTci IA=16807,IM=2147483647,IQ=127773,IR=2836,NTAB=32;
  MTci NDIV=(1+(IM-1)/NTAB);
  MTcd EPS=3.0e-16,AM=1.0/IM,RNMX=(1.0-EPS);
  static int iy=0;
  static vector<int> iv(NTAB);
  int j,k;
  double temp;

  if(idum<=0 || !iy){
    if(-idum<1) idum=1;
    else idum = -idum;
    for(j=NTAB+7;j>=0;j--){
      k=idum/IQ;
      idum=IA*(idum-k*IQ)-IR*k;
      if(idum<0) idum += IM;
      if(j<NTAB) iv[j] = idum;
    }
  iy=iv[0];
  }
  k=idum/IQ;
  idum=IA*(idum-k*IQ)-IR*k;
  if(idum<0) idum += IM;
  j=iy/NDIV;
  iy=iv[j];
  iv[j] = idum;
  if((temp=AM*iy)>RNMX) return RNMX;
  else return temp;
}



#endif /* GAUSS_H */
